Load in neccesary data

#sample_info <- read.table("sample_data.txt")
sample_info_tab <- read.table("/Users/magdalenapolder/Documents/examensarbete/scripting/sample_info.tsv", header=T, row.names=1, check.names=F, sep="\t")[-28, ]
sample_info_tab$color <- as.character(sample_info_tab$color)
count_data <- read.table("ASVs_counts.tsv", header=T, row.names=1,
             check.names=F, sep="\t")

Load neccesary packages

library(tidyverse)
library(phyloseq)
library(vegan)

Function

#Define Input
count_input <- count_data
info_input <- sample_info_tab
rarefaction_threshhold <- 50
repeat_amount <- 50


#Set up working files
rarified_count <- rrarefy(t(count_input),rarefaction_threshhold)
duplicated_info <- info_input

#Perform repeated rarefaction
if (repeat_amount > 1) {
  for (x in 2:repeat_amount){
    rarified_count <- rbind(rarified_count,rrarefy(t(count_input),2000))
    duplicated_info <- rbind(duplicated_info, info_input)
    }
  }



rare_count_phy_repeat <- otu_table(t(rarified_count), taxa_are_rows=T)

sample_info_tab_phy_repeat <- sample_data(duplicated_info)
rare_physeq_repeat <- phyloseq(rare_count_phy_repeat,sample_info_tab_phy_repeat)

vst_pcoa_repeat <- ordinate(rare_physeq_repeat, method="PCoA", distance="bray")

plot_ordination(rare_physeq_repeat, vst_pcoa_repeat,  color = "location")


coordinates <- vst_pcoa_repeat$vectors[,c(1,2)]
sum_test <- as.data.frame(rowsum(coordinates, row.names(coordinates)))
scale_factor <- length(coordinates)/length(t(sum_test))
sum_test <- as.data.frame(sum_test/scale_factor)

my_plot <- (plot_ordination(rare_physeq_repeat, vst_pcoa_repeat,justDF = TRUE))

mean_df <- my_plot[1:length(sum_test$Axis.1),]
mean_df[,1] <- sum_test$Axis.1
mean_df[,2] <- sum_test$Axis.2

#plot(x=sum_test[,1], y=sum_test[,2])

ggplot(data = my_plot, aes(x=my_plot[,1], y=my_plot[,2],color=my_plot$location)) + labs(colour = "Location", x = "NMDS1", y ="NMDS2") + geom_point(na.rm=TRUE, shape=NA) +stat_ellipse(linetype = 1,lwd = 0.8, aes(color=my_plot$location, group=my_plot$sample_id))+geom_point(data=mean_df, mapping =aes(x=Axis.1, y=Axis.2, alpha=0, color=mean_df$location)) + geom_point(alpha = 0, na.rm=TRUE)


#plot_ordination(rare_physeq_repeat, vst_pcoa_repeat,) + geom_point(shape="cross")+
#  geom_point(data=sum_test, mapping =aes(x=MDS1, y=MDS2)) +
#  stat_ellipse(linetype = 1,lwd = 0.8) #+
  #geom_point(size=1) + labs(col="location") + 
    #geom_text(aes(label=rownames(duplicated_info), hjust=0.3, vjust=-0.4), size=3)
repeat_raref <- function(count, info, repeats, threshold, colorb, shapeb) {
  #Set up working files
  rarified_count <- rrarefy(t(count),threshold)
  duplicated_info <- info
  
  #Perform repeated rarefaction
  if (repeats > 1) {
    for (x in 2:repeats){
      rarified_count <- rbind(rarified_count,rrarefy(t(count_input),2000))
      duplicated_info <- rbind(duplicated_info, info)
      }
    }

  #Convert the input into physeq objects
  rare_count_phy_repeat <- otu_table(t(rarified_count), taxa_are_rows=T)
  sample_info_tab_phy_repeat <- sample_data(duplicated_info)
  rare_physeq_repeat <- phyloseq(rare_count_phy_repeat,sample_info_tab_phy_repeat)
  
  #Perform the Ordination calculation
  vst_pcoa <- ordinate(rare_physeq_repeat, method="NMDS", distance="bray")
  
  #Create matrix with mean position for each sample
  coordinates <- vst_pcoa$points
  sum_test <- as.data.frame(rowsum(coordinates, row.names(coordinates)))
  scale_factor <- length(coordinates)/length(t(sum_test))
  sum_test <- as.data.frame(sum_test/scale_factor)
  
  #Convert ordination result into a data fram object
  my_plot <- (plot_ordination(rare_physeq_repeat, vst_pcoa,justDF = TRUE))
  
  #Add info to mean location data frame
  mean_df <- my_plot[1:length(sum_test$MDS1),]
  mean_df[,1] <- sum_test$MDS1
  mean_df[,2] <- sum_test$MDS2

  #Create the plot and print it
  finished_plot <- ggplot(data = my_plot, aes(x=my_plot[,1], y=my_plot[,2],color=my_plot$location)) + 
    labs(colour = "Location", x = "NMDS1", y ="NMDS2") + geom_point(na.rm=TRUE, shape=NA) +
    stat_ellipse(linetype = 1,lwd = 0.8, aes(color=my_plot$location, group=my_plot$sample_id)) +
    geom_point(data=mean_df, mapping =aes(x=NMDS1, y=NMDS2, alpha=0, color=mean_df$location)) +
    geom_point(alpha = 0, na.rm=TRUE)
  
  print(finished_plot)
  return(vst_pcoa)
}
complete_ordination <- repeat_raref(count_data, sample_info_tab, repeats = 20, threshold = 50, "sample_id", "location")
Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': 'x' must be numeric
repeat_list <- c(10, 15, 20, 25, 30, 35, 40, 45, 50, 55, 60, 65, 70, 80, 90, 100, 125, 150, 175, 200)
time_data <- matrix(ncol = 4, nrow = 0)
colnames(time_data) <- c("Repeat Amount", "Threshold", "Time in sec", "Time in min" )


for (x in repeat_list) {
  print(paste("Running with ",x," repeats"))
  time1 <- Sys.time()
  complete_ordination <- repeat_raref(count_data, sample_info_tab, repeats = x, threshold = 1000, "sample_id", "location")
  time2 <-Sys.time()
  time_taken_secs <- difftime(time2, time1, units="secs")
  time_taken_min <- difftime(time2, time1, units="mins")
  time_data <- rbind(time_data, c(x, 1000, time_taken_secs, time_taken_min))
}
[1] "Running with  10  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2068946 
Run 1 stress 0.2663137 
Run 2 stress 0.208834 
Run 3 stress 0.2108935 
Run 4 stress 0.2484204 
Run 5 stress 0.2472172 
Run 6 stress 0.2341156 
Run 7 stress 0.2587074 
Run 8 stress 0.2148219 
Run 9 stress 0.2122492 
Run 10 stress 0.2671656 
Run 11 stress 0.2366999 
Run 12 stress 0.2394257 
Run 13 stress 0.2363499 
Run 14 stress 0.2392263 
Run 15 stress 0.2484121 
Run 16 stress 0.2092305 
Run 17 stress 0.2329405 
Run 18 stress 0.237903 
Run 19 stress 0.2287428 
Run 20 stress 0.2107879 
*** Best solution was not repeated -- monoMDS stopping criteria:
     3: no. of iterations >= maxit
    12: stress ratio > sratmax
     5: scale factor of the gradient < sfgrmin
[1] "Running with  15  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2037621 
Run 1 stress 0.2144868 
Run 2 stress 0.2339792 
Run 3 stress 0.2076302 
Run 4 stress 0.266159 
Run 5 stress 0.2037621 
... Procrustes: rmse 1.353935e-05  max resid 8.49778e-05 
... Similar to previous best
Run 6 stress 0.2620269 
Run 7 stress 0.2037621 
... Procrustes: rmse 1.964777e-05  max resid 0.0001768671 
... Similar to previous best
Run 8 stress 0.2278356 
Run 9 stress 0.2531525 
Run 10 stress 0.2219272 
Run 11 stress 0.2404791 
Run 12 stress 0.2309629 
Run 13 stress 0.2050069 
Run 14 stress 0.2134233 
Run 15 stress 0.2485466 
Run 16 stress 0.2721227 
Run 17 stress 0.2063026 
Run 18 stress 0.2062653 
Run 19 stress 0.2663672 
Run 20 stress 0.2544041 
*** Best solution repeated 2 times
[1] "Running with  20  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2172838 
Run 1 stress 0.247867 
Run 2 stress 0.2323725 
Run 3 stress 0.215454 
... New best solution
... Procrustes: rmse 0.008543489  max resid 0.07234141 
Run 4 stress 0.2104664 
... New best solution
... Procrustes: rmse 0.009625521  max resid 0.06633488 
Run 5 stress 0.2324129 
Run 6 stress 0.2582749 
Run 7 stress 0.2092665 
... New best solution
... Procrustes: rmse 0.00585217  max resid 0.07313217 
Run 8 stress 0.2102615 
Run 9 stress 0.2482576 
Run 10 stress 0.2107992 
Run 11 stress 0.2062045 
... New best solution
... Procrustes: rmse 0.005213015  max resid 0.05727515 
Run 12 stress 0.24784 
Run 13 stress 0.2556115 
Run 14 stress 0.2319646 
Run 15 stress 0.2621871 
Run 16 stress 0.2278621 
Run 17 stress 0.2113299 
Run 18 stress 0.2736252 
Run 19 stress 0.2080111 
Run 20 stress 0.2081056 
*** Best solution was not repeated -- monoMDS stopping criteria:
     3: no. of iterations >= maxit
     4: stress ratio > sratmax
    13: scale factor of the gradient < sfgrmin
[1] "Running with  25  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2159204 
Run 1 stress 0.2109334 
... New best solution
... Procrustes: rmse 0.01345458  max resid 0.0647212 
Run 2 stress 0.2323745 
Run 3 stress 0.2096228 
... New best solution
... Procrustes: rmse 0.00694352  max resid 0.05694282 
Run 4 stress 0.2302605 
Run 5 stress 0.2333996 
Run 6 stress 0.2242031 
Run 7 stress 0.2061253 
... New best solution
... Procrustes: rmse 0.00675013  max resid 0.07979637 
Run 8 stress 0.2481123 
Run 9 stress 0.2386248 
Run 10 stress 0.2131409 
Run 11 stress 0.2100337 
Run 12 stress 0.2562804 
Run 13 stress 0.2039217 
... New best solution
... Procrustes: rmse 0.004461687  max resid 0.03551654 
Run 14 stress 0.2133143 
Run 15 stress 0.2050656 
Run 16 stress 0.2322626 
Run 17 stress 0.2455612 
Run 18 stress 0.2270098 
Run 19 stress 0.2324001 
Run 20 stress 0.2292402 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
    18: scale factor of the gradient < sfgrmin
[1] "Running with  30  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2048868 
Run 1 stress 0.2066875 
Run 2 stress 0.2502874 
Run 3 stress 0.2077659 
Run 4 stress 0.2086475 
Run 5 stress 0.227905 
Run 6 stress 0.2256868 
Run 7 stress 0.2322326 
Run 8 stress 0.2456374 
Run 9 stress 0.2153393 
Run 10 stress 0.235093 
Run 11 stress 0.2416219 
Run 12 stress 0.2322977 
Run 13 stress 0.2491814 
Run 14 stress 0.2123347 
Run 15 stress 0.2081935 
Run 16 stress 0.2073044 
Run 17 stress 0.2079442 
Run 18 stress 0.2427943 
Run 19 stress 0.2475072 
Run 20 stress 0.2485183 
*** Best solution was not repeated -- monoMDS stopping criteria:
     3: no. of iterations >= maxit
     2: stress ratio > sratmax
    15: scale factor of the gradient < sfgrmin
[1] "Running with  35  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2159627 
Run 1 stress 0.2440572 
Run 2 stress 0.2325243 
Run 3 stress 0.2165349 
Run 4 stress 0.2365486 
Run 5 stress 0.2434068 
Run 6 stress 0.2120876 
... New best solution
... Procrustes: rmse 0.01209213  max resid 0.05763703 
Run 7 stress 0.2403782 
Run 8 stress 0.2415735 
Run 9 stress 0.2331229 
Run 10 stress 0.2611486 
Run 11 stress 0.2381923 
Run 12 stress 0.215076 
Run 13 stress 0.2491559 
Run 14 stress 0.2372107 
Run 15 stress 0.2097071 
... New best solution
... Procrustes: rmse 0.005385369  max resid 0.0493466 
Run 16 stress 0.2184571 
Run 17 stress 0.2646963 
Run 18 stress 0.2270707 
Run 19 stress 0.2363402 
Run 20 stress 0.2656087 
*** Best solution was not repeated -- monoMDS stopping criteria:
     5: no. of iterations >= maxit
     1: stress ratio > sratmax
    14: scale factor of the gradient < sfgrmin
[1] "Running with  40  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2148911 
Run 1 stress 0.2341205 
Run 2 stress 0.260825 
Run 3 stress 0.2493988 
Run 4 stress 0.2370208 
Run 5 stress 0.2376256 
Run 6 stress 0.2613994 
Run 7 stress 0.2419957 
Run 8 stress 0.2400575 
Run 9 stress 0.2500377 
Run 10 stress 0.2500017 
Run 11 stress 0.248826 
Run 12 stress 0.2099825 
... New best solution
... Procrustes: rmse 0.01046906  max resid 0.05603305 
Run 13 stress 0.2580906 
Run 14 stress 0.2055665 
... New best solution
... Procrustes: rmse 0.004565162  max resid 0.04443016 
Run 15 stress 0.246378 
Run 16 stress 0.2048284 
... New best solution
... Procrustes: rmse 0.003135009  max resid 0.04922089 
Run 17 stress 0.2309908 
Run 18 stress 0.2343338 
Run 19 stress 0.2085154 
Run 20 stress 0.2419515 
*** Best solution was not repeated -- monoMDS stopping criteria:
     7: no. of iterations >= maxit
    13: scale factor of the gradient < sfgrmin
[1] "Running with  45  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2155541 
Run 1 stress 0.2456793 
Run 2 stress 0.2422816 
Run 3 stress 0.2253852 
Run 4 stress 0.2374882 
Run 5 stress 0.2448256 
Run 6 stress 0.225317 
Run 7 stress 0.25378 
Run 8 stress 0.2476409 
Run 9 stress 0.2128044 
... New best solution
... Procrustes: rmse 0.01078795  max resid 0.05583074 
Run 10 stress 0.2105192 
... New best solution
... Procrustes: rmse 0.005344977  max resid 0.04347715 
Run 11 stress 0.2049319 
... New best solution
... Procrustes: rmse 0.004572584  max resid 0.06137488 
Run 12 stress 0.2112032 
Run 13 stress 0.2445934 
Run 14 stress 0.2138173 
Run 15 stress 0.2631151 
Run 16 stress 0.2447748 
Run 17 stress 0.2468428 
Run 18 stress 0.2385953 
Run 19 stress 0.2643583 
Run 20 stress 0.205687 
*** Best solution was not repeated -- monoMDS stopping criteria:
     5: no. of iterations >= maxit
    15: scale factor of the gradient < sfgrmin
[1] "Running with  50  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2146711 
Run 1 stress 0.247094 
Run 2 stress 0.2468967 
Run 3 stress 0.225205 
Run 4 stress 0.2373683 
Run 5 stress 0.206678 
... New best solution
... Procrustes: rmse 0.008840808  max resid 0.05044393 
Run 6 stress 0.2625275 
Run 7 stress 0.2646484 
Run 8 stress 0.2435206 
Run 9 stress 0.2123092 
Run 10 stress 0.2588298 
Run 11 stress 0.2076638 
Run 12 stress 0.2057752 
... New best solution
... Procrustes: rmse 0.003290995  max resid 0.044621 
Run 13 stress 0.2624739 
Run 14 stress 0.2124649 
Run 15 stress 0.2041138 
... New best solution
... Procrustes: rmse 0.002811395  max resid 0.04392208 
Run 16 stress 0.227348 
Run 17 stress 0.2149618 
Run 18 stress 0.2481491 
Run 19 stress 0.2539749 
Run 20 stress 0.2066677 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
    18: scale factor of the gradient < sfgrmin
[1] "Running with  55  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2149181 
Run 1 stress 0.2094504 
... New best solution
... Procrustes: rmse 0.009131158  max resid 0.04576863 
Run 2 stress 0.2630586 
Run 3 stress 0.2458759 
Run 4 stress 0.2291794 
Run 5 stress 0.2083754 
... New best solution
... Procrustes: rmse 0.004126547  max resid 0.03881932 
Run 6 stress 0.231999 
Run 7 stress 0.2452261 
Run 8 stress 0.2456985 
Run 9 stress 0.2229702 
Run 10 stress 0.2455323 
Run 11 stress 0.2602408 
Run 12 stress 0.2468942 
Run 13 stress 0.4208333 
Run 14 stress 0.2473179 
Run 15 stress 0.2249847 
Run 16 stress 0.2313205 
Run 17 stress 0.2249866 
Run 18 stress 0.2061198 
... New best solution
... Procrustes: rmse 0.004316067  max resid 0.03895677 
Run 19 stress 0.2457703 
Run 20 stress 0.2578304 
*** Best solution was not repeated -- monoMDS stopping criteria:
    20: scale factor of the gradient < sfgrmin
[1] "Running with  60  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2156977 
Run 1 stress 0.2096811 
... New best solution
... Procrustes: rmse 0.007936442  max resid 0.04992578 
Run 2 stress 0.2115546 
Run 3 stress 0.4208715 
Run 4 stress 0.2078592 
... New best solution
... Procrustes: rmse 0.003649462  max resid 0.03922537 
Run 5 stress 0.2534319 
Run 6 stress 0.2261133 
Run 7 stress 0.2467298 
Run 8 stress 0.2596728 
Run 9 stress 0.2230132 
Run 10 stress 0.4208726 
Run 11 stress 0.2081206 
... Procrustes: rmse 0.00344311  max resid 0.03870427 
Run 12 stress 0.2118499 
Run 13 stress 0.2595362 
Run 14 stress 0.2571937 
Run 15 stress 0.2286469 
Run 16 stress 0.2403542 
Run 17 stress 0.2302645 
Run 18 stress 0.2134467 
Run 19 stress 0.2645 
Run 20 stress 0.2048229 
... New best solution
... Procrustes: rmse 0.002946028  max resid 0.03562683 
*** Best solution was not repeated -- monoMDS stopping criteria:
     8: no. of iterations >= maxit
    12: scale factor of the gradient < sfgrmin
[1] "Running with  65  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.215196 
Run 1 stress 0.2367409 
Run 2 stress 0.2367428 
Run 3 stress 0.2102159 
... New best solution
... Procrustes: rmse 0.008063984  max resid 0.04271201 
Run 4 stress 0.2454367 
Run 5 stress 0.205272 
... New best solution
... Procrustes: rmse 0.003309468  max resid 0.0381221 
Run 6 stress 0.2078369 
Run 7 stress 0.2299384 
Run 8 stress 0.2068666 
Run 9 stress 0.2042599 
... New best solution
... Procrustes: rmse 0.001529567  max resid 0.03292068 
Run 10 stress 0.2329084 
Run 11 stress 0.21946 
Run 12 stress 0.2469003 
Run 13 stress 0.2446621 
Run 14 stress 0.2043244 
... Procrustes: rmse 0.0004486976  max resid 0.01891409 
Run 15 stress 0.2090608 
Run 16 stress 0.2450009 
Run 17 stress 0.2489648 
Run 18 stress 0.2289682 
Run 19 stress 0.2244397 
Run 20 stress 0.2480431 
*** Best solution was not repeated -- monoMDS stopping criteria:
     5: no. of iterations >= maxit
    15: scale factor of the gradient < sfgrmin
[1] "Running with  70  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2156611 
Run 1 stress 0.2463472 
Run 2 stress 0.2547421 
Run 3 stress 0.4209413 
Run 4 stress 0.2236026 
Run 5 stress 0.2281274 
Run 6 stress 0.2575193 
Run 7 stress 0.2308123 
Run 8 stress 0.2100399 
... New best solution
... Procrustes: rmse 0.007666282  max resid 0.04136695 
Run 9 stress 0.2082502 
... New best solution
... Procrustes: rmse 0.003740293  max resid 0.03859532 
Run 10 stress 0.2437089 
Run 11 stress 0.2635594 
Run 12 stress 0.2112224 
Run 13 stress 0.2447727 
Run 14 stress 0.2377362 
Run 15 stress 0.2593871 
Run 16 stress 0.2101711 
Run 17 stress 0.2107413 
Run 18 stress 0.2121219 
Run 19 stress 0.242272 
Run 20 stress 0.2107429 
*** Best solution was not repeated -- monoMDS stopping criteria:
     4: no. of iterations >= maxit
    16: scale factor of the gradient < sfgrmin
[1] "Running with  80  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.215822 
Run 1 stress 0.2071136 
... New best solution
... Procrustes: rmse 0.006930798  max resid 0.03995484 
Run 2 stress 0.2069358 
... New best solution
... Procrustes: rmse 0.002463907  max resid 0.03575849 
Run 3 stress 0.2463148 
Run 4 stress 0.2082768 
Run 5 stress 0.2606677 
Run 6 stress 0.2110424 
Run 7 stress 0.2284578 
Run 8 stress 0.2593379 
Run 9 stress 0.258091 
Run 10 stress 0.2426367 
Run 11 stress 0.2084672 
Run 12 stress 0.263565 
Run 13 stress 0.2610747 
Run 14 stress 0.2444638 
Run 15 stress 0.2597474 
Run 16 stress 0.2257589 
Run 17 stress 0.2482234 
Run 18 stress 0.2488592 
Run 19 stress 0.2470879 
Run 20 stress 0.2466732 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
    18: scale factor of the gradient < sfgrmin
[1] "Running with  90  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2035858 
Run 1 stress 0.2376481 
Run 2 stress 0.2328945 
Run 3 stress 0.2116773 
Run 4 stress 0.2071597 
Run 5 stress 0.2318099 
Run 6 stress 0.2320072 
Run 7 stress 0.421032 
Run 8 stress 0.2457167 
Run 9 stress 0.2463062 
Run 10 stress 0.2441054 
Run 11 stress 0.2226625 
Run 12 stress 0.421033 
Run 13 stress 0.2628761 
Run 14 stress 0.2469244 
Run 15 stress 0.2448891 
Run 16 stress 0.2503415 
Run 17 stress 0.2460393 
Run 18 stress 0.421027 
Run 19 stress 0.245976 
Run 20 stress 0.2427347 
*** Best solution was not repeated -- monoMDS stopping criteria:
     2: no. of iterations >= maxit
    18: scale factor of the gradient < sfgrmin
[1] "Running with  100  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2044245 
Run 1 stress 0.2370089 
Run 2 stress 0.2327077 
Run 3 stress 0.2119449 
Run 4 stress 0.2439213 
Run 5 stress 0.2341286 
Run 6 stress 0.2093629 
Run 7 stress 0.208026 
Run 8 stress 0.2132467 
Run 9 stress 0.2102858 
Run 10 stress 0.2644138 
Run 11 stress 0.2141068 
Run 12 stress 0.2472585 
Run 13 stress 0.240561 
Run 14 stress 0.2248437 
Run 15 stress 0.2332474 
Run 16 stress 0.2604631 
Run 17 stress 0.2187768 
Run 18 stress 0.2247018 
Run 19 stress 0.2333012 
Run 20 stress 0.2134549 
*** Best solution was not repeated -- monoMDS stopping criteria:
     7: no. of iterations >= maxit
    13: scale factor of the gradient < sfgrmin
[1] "Running with  125  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2159332 
Run 1 stress 0.246082 
Run 2 stress 0.2480258 
Run 3 stress 0.2470711 
Run 4 stress 0.2294322 
Run 5 stress 0.2391229 
Run 6 stress 0.23303 
Run 7 stress 0.204853 
... New best solution
... Procrustes: rmse 0.005227148  max resid 0.03211169 
Run 8 stress 0.4211179 
Run 9 stress 0.2302452 
Run 10 stress 0.2444754 
Run 11 stress 0.2451402 
Run 12 stress 0.2119602 
Run 13 stress 0.2474627 
Run 14 stress 0.2281534 
Run 15 stress 0.2475207 
Run 16 stress 0.232335 
Run 17 stress 0.2647527 
Run 18 stress 0.2402844 
Run 19 stress 0.2322469 
Run 20 stress 0.2285175 
*** Best solution was not repeated -- monoMDS stopping criteria:
     3: no. of iterations >= maxit
    17: scale factor of the gradient < sfgrmin
[1] "Running with  150  repeats"
Square root transformation
Wisconsin double standardization
Run 0 stress 0.2157534 
Run 1 stress 0.246089 
Run 2 stress 0.2475034 
Run 3 stress 0.4211591 
Run 4 stress 0.2330081 
Run 5 stress 0.2360439 
Run 6 stress 0.2357407 
Run 7 stress 0.2268399 
Run 8 stress 0.2105982 
... New best solution
... Procrustes: rmse 0.005350735  max resid 0.03034189 
Run 9 stress 0.2587474 
Run 10 stress 0.2255851 
Run 11 stress 0.2641206 
Run 12 stress 0.2442011 
Run 13 stress 0.226946 
Run 14 stress 0.2382041 
Run 15 stress 0.4211571 
Run 16 stress 0.2261881 
Run 17 stress 0.2331532 
Run 18 stress 0.4211552 
Run 19 stress 0.2447606 
Run 20 stress 0.263867 
*** Best solution was not repeated -- monoMDS stopping criteria:
     5: no. of iterations >= maxit
    15: scale factor of the gradient < sfgrmin
[1] "Running with  175  repeats"
Square root transformation

LS0tCnRpdGxlOiAiUmVwZWF0ZWQgUmFyZWZhY3Rpb24gRnVuY3Rpb24iCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCgpMb2FkIGluIG5lY2Nlc2FyeSBkYXRhCmBgYHtyfQojc2FtcGxlX2luZm8gPC0gcmVhZC50YWJsZSgic2FtcGxlX2RhdGEudHh0IikKc2FtcGxlX2luZm9fdGFiIDwtIHJlYWQudGFibGUoIi9Vc2Vycy9tYWdkYWxlbmFwb2xkZXIvRG9jdW1lbnRzL2V4YW1lbnNhcmJldGUvc2NyaXB0aW5nL3NhbXBsZV9pbmZvLnRzdiIsIGhlYWRlcj1ULCByb3cubmFtZXM9MSwgY2hlY2submFtZXM9Riwgc2VwPSJcdCIpWy0yOCwgXQpzYW1wbGVfaW5mb190YWIkY29sb3IgPC0gYXMuY2hhcmFjdGVyKHNhbXBsZV9pbmZvX3RhYiRjb2xvcikKY291bnRfZGF0YSA8LSByZWFkLnRhYmxlKCJBU1ZzX2NvdW50cy50c3YiLCBoZWFkZXI9VCwgcm93Lm5hbWVzPTEsCiAgICAgICAgICAgICBjaGVjay5uYW1lcz1GLCBzZXA9Ilx0IikKYGBgCgpMb2FkIG5lY2Nlc2FyeSBwYWNrYWdlcwpgYGB7cn0KbGlicmFyeSh0aWR5dmVyc2UpCmxpYnJhcnkocGh5bG9zZXEpCmxpYnJhcnkodmVnYW4pCmBgYAoKRnVuY3Rpb24KYGBge3J9CiNEZWZpbmUgSW5wdXQKY291bnRfaW5wdXQgPC0gY291bnRfZGF0YQppbmZvX2lucHV0IDwtIHNhbXBsZV9pbmZvX3RhYgpyYXJlZmFjdGlvbl90aHJlc2hob2xkIDwtIDUwCnJlcGVhdF9hbW91bnQgPC0gNTAKCgojU2V0IHVwIHdvcmtpbmcgZmlsZXMKcmFyaWZpZWRfY291bnQgPC0gcnJhcmVmeSh0KGNvdW50X2lucHV0KSxyYXJlZmFjdGlvbl90aHJlc2hob2xkKQpkdXBsaWNhdGVkX2luZm8gPC0gaW5mb19pbnB1dAoKI1BlcmZvcm0gcmVwZWF0ZWQgcmFyZWZhY3Rpb24KaWYgKHJlcGVhdF9hbW91bnQgPiAxKSB7CiAgZm9yICh4IGluIDI6cmVwZWF0X2Ftb3VudCl7CiAgICByYXJpZmllZF9jb3VudCA8LSByYmluZChyYXJpZmllZF9jb3VudCxycmFyZWZ5KHQoY291bnRfaW5wdXQpLDIwMDApKQogICAgZHVwbGljYXRlZF9pbmZvIDwtIHJiaW5kKGR1cGxpY2F0ZWRfaW5mbywgaW5mb19pbnB1dCkKICAgIH0KICB9CgoKCnJhcmVfY291bnRfcGh5X3JlcGVhdCA8LSBvdHVfdGFibGUodChyYXJpZmllZF9jb3VudCksIHRheGFfYXJlX3Jvd3M9VCkKCnNhbXBsZV9pbmZvX3RhYl9waHlfcmVwZWF0IDwtIHNhbXBsZV9kYXRhKGR1cGxpY2F0ZWRfaW5mbykKcmFyZV9waHlzZXFfcmVwZWF0IDwtIHBoeWxvc2VxKHJhcmVfY291bnRfcGh5X3JlcGVhdCxzYW1wbGVfaW5mb190YWJfcGh5X3JlcGVhdCkKCnZzdF9wY29hX3JlcGVhdCA8LSBvcmRpbmF0ZShyYXJlX3BoeXNlcV9yZXBlYXQsIG1ldGhvZD0iUENvQSIsIGRpc3RhbmNlPSJicmF5IikKCnBsb3Rfb3JkaW5hdGlvbihyYXJlX3BoeXNlcV9yZXBlYXQsIHZzdF9wY29hX3JlcGVhdCwgIGNvbG9yID0gImxvY2F0aW9uIikKCmNvb3JkaW5hdGVzIDwtIHZzdF9wY29hX3JlcGVhdCR2ZWN0b3JzWyxjKDEsMildCnN1bV90ZXN0IDwtIGFzLmRhdGEuZnJhbWUocm93c3VtKGNvb3JkaW5hdGVzLCByb3cubmFtZXMoY29vcmRpbmF0ZXMpKSkKc2NhbGVfZmFjdG9yIDwtIGxlbmd0aChjb29yZGluYXRlcykvbGVuZ3RoKHQoc3VtX3Rlc3QpKQpzdW1fdGVzdCA8LSBhcy5kYXRhLmZyYW1lKHN1bV90ZXN0L3NjYWxlX2ZhY3RvcikKCm15X3Bsb3QgPC0gKHBsb3Rfb3JkaW5hdGlvbihyYXJlX3BoeXNlcV9yZXBlYXQsIHZzdF9wY29hX3JlcGVhdCxqdXN0REYgPSBUUlVFKSkKCm1lYW5fZGYgPC0gbXlfcGxvdFsxOmxlbmd0aChzdW1fdGVzdCRBeGlzLjEpLF0KbWVhbl9kZlssMV0gPC0gc3VtX3Rlc3QkQXhpcy4xCm1lYW5fZGZbLDJdIDwtIHN1bV90ZXN0JEF4aXMuMgoKI3Bsb3QoeD1zdW1fdGVzdFssMV0sIHk9c3VtX3Rlc3RbLDJdKQoKZ2dwbG90KGRhdGEgPSBteV9wbG90LCBhZXMoeD1teV9wbG90WywxXSwgeT1teV9wbG90WywyXSxjb2xvcj1teV9wbG90JGxvY2F0aW9uKSkgKyBsYWJzKGNvbG91ciA9ICJMb2NhdGlvbiIsIHggPSAiTk1EUzEiLCB5ID0iTk1EUzIiKSArIGdlb21fcG9pbnQobmEucm09VFJVRSwgc2hhcGU9TkEpICtzdGF0X2VsbGlwc2UobGluZXR5cGUgPSAxLGx3ZCA9IDAuOCwgYWVzKGNvbG9yPW15X3Bsb3QkbG9jYXRpb24sIGdyb3VwPW15X3Bsb3Qkc2FtcGxlX2lkKSkrZ2VvbV9wb2ludChkYXRhPW1lYW5fZGYsIG1hcHBpbmcgPWFlcyh4PUF4aXMuMSwgeT1BeGlzLjIsIGFscGhhPTAsIGNvbG9yPW1lYW5fZGYkbG9jYXRpb24pKSArIGdlb21fcG9pbnQoYWxwaGEgPSAwLCBuYS5ybT1UUlVFKQoKI3Bsb3Rfb3JkaW5hdGlvbihyYXJlX3BoeXNlcV9yZXBlYXQsIHZzdF9wY29hX3JlcGVhdCwpICsgZ2VvbV9wb2ludChzaGFwZT0iY3Jvc3MiKSsKIyAgZ2VvbV9wb2ludChkYXRhPXN1bV90ZXN0LCBtYXBwaW5nID1hZXMoeD1NRFMxLCB5PU1EUzIpKSArCiMgIHN0YXRfZWxsaXBzZShsaW5ldHlwZSA9IDEsbHdkID0gMC44KSAjKwogICNnZW9tX3BvaW50KHNpemU9MSkgKyBsYWJzKGNvbD0ibG9jYXRpb24iKSArIAogICAgI2dlb21fdGV4dChhZXMobGFiZWw9cm93bmFtZXMoZHVwbGljYXRlZF9pbmZvKSwgaGp1c3Q9MC4zLCB2anVzdD0tMC40KSwgc2l6ZT0zKQoKYGBgCgpgYGB7cn0KcmVwZWF0X3JhcmVmIDwtIGZ1bmN0aW9uKGNvdW50LCBpbmZvLCByZXBlYXRzLCB0aHJlc2hvbGQsIGNvbG9yYiwgc2hhcGViKSB7CiAgI1NldCB1cCB3b3JraW5nIGZpbGVzCiAgcmFyaWZpZWRfY291bnQgPC0gcnJhcmVmeSh0KGNvdW50KSx0aHJlc2hvbGQpCiAgZHVwbGljYXRlZF9pbmZvIDwtIGluZm8KICAKICAjUGVyZm9ybSByZXBlYXRlZCByYXJlZmFjdGlvbgogIGlmIChyZXBlYXRzID4gMSkgewogICAgZm9yICh4IGluIDI6cmVwZWF0cyl7CiAgICAgIHJhcmlmaWVkX2NvdW50IDwtIHJiaW5kKHJhcmlmaWVkX2NvdW50LHJyYXJlZnkodChjb3VudF9pbnB1dCksMjAwMCkpCiAgICAgIGR1cGxpY2F0ZWRfaW5mbyA8LSByYmluZChkdXBsaWNhdGVkX2luZm8sIGluZm8pCiAgICAgIH0KICAgIH0KCiAgI0NvbnZlcnQgdGhlIGlucHV0IGludG8gcGh5c2VxIG9iamVjdHMKICByYXJlX2NvdW50X3BoeV9yZXBlYXQgPC0gb3R1X3RhYmxlKHQocmFyaWZpZWRfY291bnQpLCB0YXhhX2FyZV9yb3dzPVQpCiAgc2FtcGxlX2luZm9fdGFiX3BoeV9yZXBlYXQgPC0gc2FtcGxlX2RhdGEoZHVwbGljYXRlZF9pbmZvKQogIHJhcmVfcGh5c2VxX3JlcGVhdCA8LSBwaHlsb3NlcShyYXJlX2NvdW50X3BoeV9yZXBlYXQsc2FtcGxlX2luZm9fdGFiX3BoeV9yZXBlYXQpCiAgCiAgI1BlcmZvcm0gdGhlIE9yZGluYXRpb24gY2FsY3VsYXRpb24KICB2c3RfcGNvYSA8LSBvcmRpbmF0ZShyYXJlX3BoeXNlcV9yZXBlYXQsIG1ldGhvZD0iTk1EUyIsIGRpc3RhbmNlPSJicmF5IikKICAKICAjQ3JlYXRlIG1hdHJpeCB3aXRoIG1lYW4gcG9zaXRpb24gZm9yIGVhY2ggc2FtcGxlCiAgY29vcmRpbmF0ZXMgPC0gdnN0X3Bjb2EkcG9pbnRzCiAgc3VtX3Rlc3QgPC0gYXMuZGF0YS5mcmFtZShyb3dzdW0oY29vcmRpbmF0ZXMsIHJvdy5uYW1lcyhjb29yZGluYXRlcykpKQogIHNjYWxlX2ZhY3RvciA8LSBsZW5ndGgoY29vcmRpbmF0ZXMpL2xlbmd0aCh0KHN1bV90ZXN0KSkKICBzdW1fdGVzdCA8LSBhcy5kYXRhLmZyYW1lKHN1bV90ZXN0L3NjYWxlX2ZhY3RvcikKICAKICAjQ29udmVydCBvcmRpbmF0aW9uIHJlc3VsdCBpbnRvIGEgZGF0YSBmcmFtIG9iamVjdAogIG15X3Bsb3QgPC0gKHBsb3Rfb3JkaW5hdGlvbihyYXJlX3BoeXNlcV9yZXBlYXQsIHZzdF9wY29hLGp1c3RERiA9IFRSVUUpKQogIAogICNBZGQgaW5mbyB0byBtZWFuIGxvY2F0aW9uIGRhdGEgZnJhbWUKICBtZWFuX2RmIDwtIG15X3Bsb3RbMTpsZW5ndGgoc3VtX3Rlc3QkTURTMSksXQogIG1lYW5fZGZbLDFdIDwtIHN1bV90ZXN0JE1EUzEKICBtZWFuX2RmWywyXSA8LSBzdW1fdGVzdCRNRFMyCgogICNDcmVhdGUgdGhlIHBsb3QgYW5kIHByaW50IGl0CiAgZmluaXNoZWRfcGxvdCA8LSBnZ3Bsb3QoZGF0YSA9IG15X3Bsb3QsIGFlcyh4PW15X3Bsb3RbLDFdLCB5PW15X3Bsb3RbLDJdLGNvbG9yPW15X3Bsb3QkbG9jYXRpb24pKSArIAogICAgbGFicyhjb2xvdXIgPSAiTG9jYXRpb24iLCB4ID0gIk5NRFMxIiwgeSA9Ik5NRFMyIikgKyBnZW9tX3BvaW50KG5hLnJtPVRSVUUsIHNoYXBlPU5BKSArCiAgICBzdGF0X2VsbGlwc2UobGluZXR5cGUgPSAxLGx3ZCA9IDAuOCwgYWVzKGNvbG9yPW15X3Bsb3QkbG9jYXRpb24sIGdyb3VwPW15X3Bsb3Qkc2FtcGxlX2lkKSkgKwogICAgZ2VvbV9wb2ludChkYXRhPW1lYW5fZGYsIG1hcHBpbmcgPWFlcyh4PU5NRFMxLCB5PU5NRFMyLCBhbHBoYT0wLCBjb2xvcj1tZWFuX2RmJGxvY2F0aW9uKSkgKwogICAgZ2VvbV9wb2ludChhbHBoYSA9IDAsIG5hLnJtPVRSVUUpCiAgCiAgcHJpbnQoZmluaXNoZWRfcGxvdCkKICByZXR1cm4odnN0X3Bjb2EpCn0KYGBgCgpgYGB7cn0KY29tcGxldGVfb3JkaW5hdGlvbiA8LSByZXBlYXRfcmFyZWYoY291bnRfZGF0YSwgc2FtcGxlX2luZm9fdGFiLCByZXBlYXRzID0gMjAsIHRocmVzaG9sZCA9IDUwLCAic2FtcGxlX2lkIiwgImxvY2F0aW9uIikKCmBgYAoKYGBge3J9CnJlcGVhdF9saXN0IDwtIGMoMTAsIDE1LCAyMCwgMjUsIDMwLCAzNSwgNDAsIDQ1LCA1MCwgNTUsIDYwLCA2NSwgNzAsIDgwLCA5MCwgMTAwLCAxMjUsIDE1MCwgMTc1LCAyMDApCnRpbWVfZGF0YSA8LSBtYXRyaXgobmNvbCA9IDQsIG5yb3cgPSAwKQpjb2xuYW1lcyh0aW1lX2RhdGEpIDwtIGMoIlJlcGVhdCBBbW91bnQiLCAiVGhyZXNob2xkIiwgIlRpbWUgaW4gc2VjIiwgIlRpbWUgaW4gbWluIiApCgoKZm9yICh4IGluIHJlcGVhdF9saXN0KSB7CiAgcHJpbnQocGFzdGUoIlJ1bm5pbmcgd2l0aCAiLHgsIiByZXBlYXRzIikpCiAgdGltZTEgPC0gU3lzLnRpbWUoKQogIGNvbXBsZXRlX29yZGluYXRpb24gPC0gcmVwZWF0X3JhcmVmKGNvdW50X2RhdGEsIHNhbXBsZV9pbmZvX3RhYiwgcmVwZWF0cyA9IHgsIHRocmVzaG9sZCA9IDEwMDAsICJzYW1wbGVfaWQiLCAibG9jYXRpb24iKQogIHRpbWUyIDwtU3lzLnRpbWUoKQogIHRpbWVfdGFrZW5fc2VjcyA8LSBkaWZmdGltZSh0aW1lMiwgdGltZTEsIHVuaXRzPSJzZWNzIikKICB0aW1lX3Rha2VuX21pbiA8LSBkaWZmdGltZSh0aW1lMiwgdGltZTEsIHVuaXRzPSJtaW5zIikKICB0aW1lX2RhdGEgPC0gcmJpbmQodGltZV9kYXRhLCBjKHgsIDEwMDAsIHRpbWVfdGFrZW5fc2VjcywgdGltZV90YWtlbl9taW4pKQp9CgoKYGBgCgpgYGB7cn0KcGxvdCh4ID0gdGltZV9kYXRhWywxXSwgeSA9IHRpbWVfZGF0YVssNF0pCgpjb2xuYW1lcyh0aW1lX2RhdGEpIDwtIGMoIlJlcGVhdF9BbW91bnQiLCAiVGhyZXNob2xkIiwgIlRpbWUgaW4gc2VjIiwgIlRpbWVfaW5fbWluIiApCnRpbWVfZGF0YSA8LSBhcy5kYXRhLmZyYW1lKHRpbWVfZGF0YSkKCmdncGxvdCh0aW1lX2RhdGEsYWVzKHg9UmVwZWF0X0Ftb3VudCx5PVRpbWVfaW5fbWluKSkgKyBnZW9tX3BvaW50KCkgKyB4bGltKDAsMTAwMCkgKwpzdGF0X3Ntb290aChtZXRob2Q9ImdhbSIsZnVsbHJhbmdlPVRSVUUpCgoKZ2dwbG90KHRpbWVfZGF0YSxhZXMoeD1SZXBlYXRfQW1vdW50LHk9VGltZV9pbl9taW4pKSsKICBnZW9tX3BvaW50KCkrCiAgc3RhdF9zbW9vdGgobWV0aG9kID0gImxtIiwgZm9ybXVsYSA9IHkgfiB4ICsgSSh4XjIpLCBzaXplID0gMSwgZnVsbHJhbmdlID0gVCkreGxpbSgwLDEwMDApCmBgYAoK